library(tidyverse)

library(broom)     # tidy model data
library(patchwork) # combine ggplots
library(ggthemes)  # ggplot Stata theme
library(estimatr)  # robust standard-errors
##  tidy lm model results with broom (and ggplot)

tidy_lm <- function(fit) {
  tidy(fit) %>% mutate_if(is.numeric, round, 2)
}

glance_lm <- function(fit) {
  glance(fit) %>% 
    select(r.squared, adj.r.squared, statistic, p.value, df, nobs) %>% 
    mutate_if(is.numeric, round, 2)
}

plot_lm <- function(fit,   bw = 25) {
  dt_pl<-  
    augment(fit) %>% 
    mutate(einkommen_geschaetzt = .fitted,
           residuen = .resid)

  pl1 <- 
    ggplot(dt_pl, aes(x = residuen)) +
    geom_histogram(binwidth = bw, colour = "black", fill = "grey") +
    stat_function(
      fun = function(x) dnorm(x, mean = mean(dt_pl$residuen),
                              sd = sd(dt_pl$residuen))
                        * bw * nrow(dt_pl),
      colour = "blue"
    ) +
   theme_stata()

  pl2 <- 
    ggplot(data = dt_pl, aes(x = einkommen_geschaetzt, y = residuen)) +
    geom_point(alpha = 0.2) +
    geom_smooth(se = FALSE) +
    theme_stata()

  pl1 + pl2 + plot_layout(widths = c(2, 3))
}

Wunsch I · Perfekte Daten

In der Vorlesung 🧑‍🏫 am 19. November 2020 gab es eine Fee 🧚 die Wünsche 🧞 erfüllt.

Den Wunsch 🪄 nach perfekten Daten 📈 erfüllt sie nicht 😥.

Die Fee 🧚 überlegt, nutzt R 👩‍💻 und schafft fast perfekte Daten 📈.

param_const <- 500
param_alter <- 25

n <- 1000
sd_einkommen <-  100

dt_raw <- 
  tibble( 
    id = 1:n,
    alter = seq(20, 60, length.out = n) %>% round(2),
    einkommen_linear = round(param_const + param_alter*alter)
    ) %>% 
  mutate(einkommen = round(einkommen_linear + rnorm(n, sd = sd_einkommen)))

\[ Einkommen_{ID} = 500 + 25 \cdot Alter_{ID} \: + ( \: Fehler_{ID} \: )\]

ggplot(dt_raw, aes(x = alter, y = einkommen)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = lm) +
  theme_stata()

Wunsch II · Diagnose für SoWis

Die Studierenden 🙇🙇🙇🙇 hatten einen eigenen Wunsch 🪄 an die Fee 🧚.

Kannst Du uns Regressions-Diagnose 📈 mit einer sozialwissenschaftlichen Geschichte 📚 erklären?

Daten Fee

Wir stellen uns Postbeamte 📯 vor die regelmäßig Lohnerhöhungen 💰 strikt nach dem Senioritätsprinzip 🥸 erhalten.

Wir haben also ein Modell 📈 bei dem das Einkommen 💰 vom Alter 🧑👴 abhängt.

Das wahre Einkommen 💰 (einkommen_linear) ist uns bekannt und unsere beobachteten Daten 🕵️ (einkommen) weichen davon normal-verteilt 🧮 ab

So hat die Fee 🧚 uns die Daten 📊 gegeben.

ggplot(dt_raw, aes(x = alter, y = einkommen)) +
  geom_point(alpha = 0.2) +
  geom_jitter(aes(y = einkommen_linear),
              width = 0.1, alpha = 0.5,
              colour = "darkgreen", shape = ".") +
  theme_stata()


DT::datatable(dt_raw)

Modell (fast) perfekt

Wir bekommen schöne Ergebnisse bei der Schätzung des Modells. 😄

fit <- lm(einkommen ~ alter, data = dt_raw)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit)

Annahmen (1-7)

  1. korrekte Spezifikation (Linearität)
  2. Erwartungswert der Residuen (lokal) = 0
  3. Normalverteilte Residuen
  4. Varianzhomogenität (gleichmäßige Streuung Residuen)
  5. Berücksichtigung einflussreicher Fälle (keine Ausreisser)
  6. keine Multikollinearität (geringe Korrelation unabhängiger Variablen)
  7. keine Autokorrelation der Residuen (insbes. Zeitreihen)

Korrekte Spezifikation

Dritt-Variablen

Briefträger ✉️ arbeiten oft nur halbtags 🕐 und verdienen 💰 dann auch nur die Hälfte.

dt <-
  dt_raw %>% 
  mutate(einkommen = if_else(id %% 3 == 0, einkommen/2 + rnorm(n, sd = 50), einkommen),
         halbtags = if_else(id %% 3 == 0, "ja", "nein") %>% fct_relevel("ja", after = Inf))

ggplot(dt, aes(x = alter, y = einkommen, colour = halbtags)) +
  geom_point(alpha = 0.2) +
  theme_stata()

Ohne Dritt-Variable 😟

fit <- lm(einkommen ~ alter, data = dt)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit, bw=50)

Mit Dritt-Variable 🙂

fit <- lm(einkommen ~ alter + halbtags, data = dt)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit)

Noch besser mit Dritt-Variable und Interaktions-Effekt 😃

fit <- lm(einkommen ~ alter*halbtags, data = dt)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit)

Funktionale Form

Bei der Bundespost 📯 wurde eine Bezahlung 💰 nach Humankapital eingeführt.

Es wurde erkannt, dass bei jungen Beamten 🧑 die Produktivität stärker wächst als bei Älteren 🧓. Daher erhalten sie höhere Lohnzuwächse 💰.

dt <-
  dt_raw %>% 
  mutate(einkommen_polynom = -300 + 80*alter - 0.725*alter^2,
         einkommen = einkommen_polynom + rnorm(n, sd = sd_einkommen))

ggplot(dt, aes(x = alter, y = einkommen)) +
  geom_point(alpha = 0.2) +
    geom_jitter(aes(y = einkommen_polynom),
                width = 0.1, alpha = 0.5,
                colour = "darkblue", shape = ".") +
  theme_stata()

Ein lineares Modell ist falsch spezifiziert. 😟 Die Regressions-Diagnose zeigt dies. ☝️

fit <- lm(einkommen ~ alter, data = dt)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit)

Korrekt ist ein Modell mit einem Polynom (hier einem quadratischem Term). 😃

fit <- lm(einkommen ~ alter + I(alter^2), data = dt)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit)

( Den R Hinweis auf Multi-Kollinearität ignorieren wir einfach. Natürlich gibt es einen Zusammenhang zwischen alter und dem quadrierten alter^2. )

Varianzhomogenität

Bei der Post 📯 haben hohe Beamte und Minster ein höheres Einkommen 💰 als einfache Beamte. Aber nicht alle werden hohe Beamte 🧑‍💼 oder Post-Minister 🏤.

dt <- dt_raw %>% mutate(einkommen = einkommen_linear + id/500 * rnorm(n, sd = sd_einkommen))

ggplot(dt, aes(x = alter, y = einkommen)) +
  geom_point(alpha = 0.2) +
    geom_jitter(aes(y = einkommen_linear),
                width = 0.1, alpha = 0.5,
                colour = "darkgreen", shape = ".") +
  theme_stata()

Im Standard-Modell erhalten wir verzerrte Standard-Fehler. 😟

fit <- lm(einkommen ~ alter, data = dt)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit)

Besser ist ein Modell mit robusten 💪 Standard-Fehlern. 😄

fit <- lm_robust(einkommen ~ alter, data = dt)
tidy_lm(fit)
glance(fit) %>% mutate_if(is.numeric, round, 2) 

Multikollinearität

Ältere 👴👵 Postbeamte 📯 sind kleiner als Jüngere 👧👦.

dt <- dt_raw %>% mutate(groesse_cm = round(seq(180, 160, length.out = n) + rnorm(n, sd = 2.5), 3))

ggplot(dt, aes(x = alter, y = groesse_cm)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = lm) +
  theme_stata()

Es gibt dann auch einen Zusammenhang 📉 zwischen groesse_cm 👵 und einkommen 💰, natürlich keinen kausalen ⚙️.

ggplot(dt, aes(x = groesse_cm, y = einkommen)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = lm) +
  theme_stata()

fit <- lm(einkommen ~ groesse_cm, data = dt)
tidy_lm(fit)
glance_lm(fit)

Nehmen wir alter 👵 und groesse_cm 📐 in einem Modell 📉 auf haben wir Multikollinearität 🔗 😟

fit <- lm(einkommen ~ alter + groesse_cm, data = dt)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit)

Multikollinearität im Modell 📉 erkennen 🕵️ wir auch an der sehr starken Korrellation der unabhängigen Variablen.

dt %>% 
  select(alter, einkommen) %>% 
  cor() %>% 
  as_tibble() %>% 
  mutate_all(round, 2)

Ausreisser

Bei der Datenübertragung 📡 gab es einen Fehler ⚠️ und die letzte Ziffer 🔢 wurde bei einigen Gehaltsangaben 💰 entfernt.

id_random <- runif(80, min = 1, max = n) %>% round()

dt <-
  dt_raw %>% 
  mutate(einkommen = if_else(id %in% id_random, round(einkommen/10), einkommen))

ggplot(dt, aes(x = alter, y = einkommen)) +
  geom_point(alpha = 0.2) +
  theme_stata()

Wir erhalten eine verzerrte Schätzung. 😟

fit <- lm(einkommen ~ alter, data = dt)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit, bw=50)

Die falschen Beobachtungen müssen wir entfernen (oder korrigieren). 😃

dt_tmp <- dt %>% filter(einkommen > 500)

fit <- lm(einkommen ~ alter, data = dt_tmp)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit)

Autokorrelation

Die Postbeamten 📯 wurden 2010 und 2020 📅 nach ihrem Einkommen 💰 befragt. ☎️

dt_tmp <-
  dt_raw %>% 
  mutate(einkommen = round(einkommen - param_alter*10 + rnorm(n, sd = 10)),
         alter = alter - 10,
         befragung = 2010) %>% 
  filter(alter >= 20)

dt <-
  dt_raw %>% 
  mutate(befragung = 2020) %>% 
  rbind(dt_tmp) %>% 
  arrange(id, befragung) %>% 
  relocate(befragung, .after = id)

DT::datatable(dt %>% filter(id > 500))  # Daten ab ID 500 um Panel-Struktur zu verdeutlichen 

Postbeamte 📯 die 2010 📅 älter als 50 waren sind inzwischen in Pension 👵 und haben an der 2020 Befragung ☎️ nicht mehr teilgenommen.

dt_tmp <- dt %>% mutate(befragung = factor(befragung))

ggplot(dt_tmp, aes(x = alter, y = einkommen, colour = befragung)) +
  geom_point(alpha = 0.2) +
  theme_stata()

Im Modell (n = 1750) sind Beobachtungen (und Residuen) nicht mehr unabhängig. 😟 Wir haben ja viele zweimal befragt. ☎️

fit <- lm(einkommen ~ alter, data = dt_tmp)
tidy_lm(fit)
glance_lm(fit)
plot_lm(fit)

Aufgabe

Diskutieren Sie 🙇🙇🙇 die Regressions-Annahmen 📈 und zeigen Sie ob und wie Verletzungen 🏥 der Annahmen erkannt 🕵️ werden können.


LS0tCnRpdGxlOiAiUmVncmVzc2lvbnMtRGlhZ25vc2Ugwrcg8J+UjSIKb3V0cHV0OgogIGh0bWxfbm90ZWJvb2s6CiAgICBjb2RlX2ZvbGRpbmc6IGhpZGUKICAgIHRvYzogeWVzCiAgICB0b2NfZmxvYXQ6IHllcwotLS0KCgpgYGB7ciwgbWVzc2FnZT1GQUxTRX0KbGlicmFyeSh0aWR5dmVyc2UpCgpsaWJyYXJ5KGJyb29tKSAgICAgIyB0aWR5IG1vZGVsIGRhdGEKbGlicmFyeShwYXRjaHdvcmspICMgY29tYmluZSBnZ3Bsb3RzCmxpYnJhcnkoZ2d0aGVtZXMpICAjIGdncGxvdCBTdGF0YSB0aGVtZQpsaWJyYXJ5KGVzdGltYXRyKSAgIyByb2J1c3Qgc3RhbmRhcmQtZXJyb3JzCmBgYAoKCmBgYHtyLCBtZXNzYWdlPUZBTFNFfQojIyAgdGlkeSBsbSBtb2RlbCByZXN1bHRzIHdpdGggYnJvb20gKGFuZCBnZ3Bsb3QpCgp0aWR5X2xtIDwtIGZ1bmN0aW9uKGZpdCkgewogIHRpZHkoZml0KSAlPiUgbXV0YXRlX2lmKGlzLm51bWVyaWMsIHJvdW5kLCAyKQp9CgpnbGFuY2VfbG0gPC0gZnVuY3Rpb24oZml0KSB7CiAgZ2xhbmNlKGZpdCkgJT4lIAogICAgc2VsZWN0KHIuc3F1YXJlZCwgYWRqLnIuc3F1YXJlZCwgc3RhdGlzdGljLCBwLnZhbHVlLCBkZiwgbm9icykgJT4lIAogICAgbXV0YXRlX2lmKGlzLm51bWVyaWMsIHJvdW5kLCAyKQp9CgpwbG90X2xtIDwtIGZ1bmN0aW9uKGZpdCwgICBidyA9IDI1KSB7CiAgZHRfcGw8LSAgCiAgICBhdWdtZW50KGZpdCkgJT4lIAogICAgbXV0YXRlKGVpbmtvbW1lbl9nZXNjaGFldHp0ID0gLmZpdHRlZCwKICAgICAgICAgICByZXNpZHVlbiA9IC5yZXNpZCkKCiAgcGwxIDwtIAogICAgZ2dwbG90KGR0X3BsLCBhZXMoeCA9IHJlc2lkdWVuKSkgKwogICAgZ2VvbV9oaXN0b2dyYW0oYmlud2lkdGggPSBidywgY29sb3VyID0gImJsYWNrIiwgZmlsbCA9ICJncmV5IikgKwogICAgc3RhdF9mdW5jdGlvbigKICAgICAgZnVuID0gZnVuY3Rpb24oeCkgZG5vcm0oeCwgbWVhbiA9IG1lYW4oZHRfcGwkcmVzaWR1ZW4pLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICBzZCA9IHNkKGR0X3BsJHJlc2lkdWVuKSkKICAgICAgICAgICAgICAgICAgICAgICAgKiBidyAqIG5yb3coZHRfcGwpLAogICAgICBjb2xvdXIgPSAiYmx1ZSIKICAgICkgKwogICB0aGVtZV9zdGF0YSgpCgogIHBsMiA8LSAKICAgIGdncGxvdChkYXRhID0gZHRfcGwsIGFlcyh4ID0gZWlua29tbWVuX2dlc2NoYWV0enQsIHkgPSByZXNpZHVlbikpICsKICAgIGdlb21fcG9pbnQoYWxwaGEgPSAwLjIpICsKICAgIGdlb21fc21vb3RoKHNlID0gRkFMU0UpICsKICAgIHRoZW1lX3N0YXRhKCkKCiAgcGwxICsgcGwyICsgcGxvdF9sYXlvdXQod2lkdGhzID0gYygyLCAzKSkKfQpgYGAKCgojIFd1bnNjaCBJIMK3IFBlcmZla3RlIERhdGVuCgpJbiBkZXIgVm9ybGVzdW5nIPCfp5HigI3wn4+rIGFtIDE5LiBOb3ZlbWJlciAyMDIwIGdhYiBlcyBlaW5lIEZlZSDwn6eaIGRpZSBXw7xuc2NoZSDwn6eeIGVyZsO8bGx0LgoKIVtdKGJpbGRlci9mYWlyLmpwZykKCkRlbiBXdW5zY2gg8J+qhCBuYWNoIHBlcmZla3RlbiBEYXRlbiDwn5OIIGVyZsO8bGx0IHNpZSBuaWNodCDwn5ilLgogCkRpZSBGZWUg8J+nmiDDvGJlcmxlZ3QsIG51dHp0IFtfX1JfX10oaHR0cHM6Ly9lZHVjYXRpb24ucnN0dWRpby5jb20vbGVhcm4vYmVnaW5uZXIvKSDwn5Gp4oCN8J+SuyB1bmQgc2NoYWZmdCBmYXN0IHBlcmZla3RlIERhdGVuIPCfk4guCgpgYGB7cn0KcGFyYW1fY29uc3QgPC0gNTAwCnBhcmFtX2FsdGVyIDwtIDI1CgpuIDwtIDEwMDAKc2RfZWlua29tbWVuIDwtICAxMDAKCmR0X3JhdyA8LSAKICB0aWJibGUoIAogICAgaWQgPSAxOm4sCiAgICBhbHRlciA9IHNlcSgyMCwgNjAsIGxlbmd0aC5vdXQgPSBuKSAlPiUgcm91bmQoMiksCiAgICBlaW5rb21tZW5fbGluZWFyID0gcm91bmQocGFyYW1fY29uc3QgKyBwYXJhbV9hbHRlciphbHRlcikKICAgICkgJT4lIAogIG11dGF0ZShlaW5rb21tZW4gPSByb3VuZChlaW5rb21tZW5fbGluZWFyICsgcm5vcm0obiwgc2QgPSBzZF9laW5rb21tZW4pKSkKYGBgCgoKJCQgRWlua29tbWVuX3tJRH0gPSA1MDAgKyAyNSBcY2RvdCBBbHRlcl97SUR9IFw6ICsgKCBcOiBGZWhsZXJfe0lEfSBcOiApJCQKCmBgYHtyfQpnZ3Bsb3QoZHRfcmF3LCBhZXMoeCA9IGFsdGVyLCB5ID0gZWlua29tbWVuKSkgKwogIGdlb21fcG9pbnQoYWxwaGEgPSAwLjIpICsKICBnZW9tX3Ntb290aChtZXRob2QgPSBsbSkgKwogIHRoZW1lX3N0YXRhKCkKYGBgCgoKIyBXdW5zY2ggSUkgwrcgRGlhZ25vc2UgZsO8ciBTb1dpcwoKRGllIFN0dWRpZXJlbmRlbiDwn5mH8J+Zh/CfmYfwn5mHIGhhdHRlbiBlaW5lbiBlaWdlbmVuIFd1bnNjaCDwn6qEIGFuIGRpZSBGZWUg8J+nmi4KCkthbm5zdCBEdSB1bnMgUmVncmVzc2lvbnMtRGlhZ25vc2Ug8J+TiCBtaXQgZWluZXIgc296aWFsd2lzc2Vuc2NoYWZ0bGljaGVuIEdlc2NoaWNodGUg8J+TmiBlcmtsw6RyZW4/CgohW10oYmlsZGVyL2ZhaXIuanBnKQoKIyMgRGF0ZW4gRmVlCgpXaXIgc3RlbGxlbiB1bnMgUG9zdGJlYW10ZSDwn5OvIHZvciBkaWUgcmVnZWxtw6TDn2lnIExvaG5lcmjDtmh1bmdlbiDwn5KwIHN0cmlrdCBuYWNoIGRlbSBTZW5pb3JpdMOkdHNwcmluemlwIPCfpbggZXJoYWx0ZW4uIAoKV2lyIGhhYmVuIGFsc28gZWluIE1vZGVsbCDwn5OIIGJlaSBkZW0gZGFzIF9fRWlua29tbWVuX18g8J+SsCB2b20gX19BbHRlcl9fIPCfp5Hwn5G0IGFiaMOkbmd0LgoKIVtdKGJpbGRlci9wb3N0YmVhbXRlLmpwZykKCkRhcyB3YWhyZSBFaW5rb21tZW4g8J+SsCAoYGVpbmtvbW1lbl9saW5lYXJgKSBpc3QgdW5zIGJla2FubnQgdW5kIHVuc2VyZSBiZW9iYWNodGV0ZW4gRGF0ZW4g8J+Vte+4jyAoYGVpbmtvbW1lbmApIHdlaWNoZW4gZGF2b24gbm9ybWFsLXZlcnRlaWx0IPCfp64gYWIKClNvIGhhdCBkaWUgRmVlIPCfp5ogdW5zIGRpZSBEYXRlbiDwn5OKIGdlZ2ViZW4uCgpgYGB7cn0KZ2dwbG90KGR0X3JhdywgYWVzKHggPSBhbHRlciwgeSA9IGVpbmtvbW1lbikpICsKICBnZW9tX3BvaW50KGFscGhhID0gMC4yKSArCiAgZ2VvbV9qaXR0ZXIoYWVzKHkgPSBlaW5rb21tZW5fbGluZWFyKSwKICAgICAgICAgICAgICB3aWR0aCA9IDAuMSwgYWxwaGEgPSAwLjUsCiAgICAgICAgICAgICAgY29sb3VyID0gImRhcmtncmVlbiIsIHNoYXBlID0gIi4iKSArCiAgdGhlbWVfc3RhdGEoKQoKRFQ6OmRhdGF0YWJsZShkdF9yYXcpCmBgYAoKIyMgTW9kZWxsIChmYXN0KSBwZXJmZWt0CgpXaXIgYmVrb21tZW4gc2Now7ZuZSBFcmdlYm5pc3NlIGJlaSBkZXIgU2Now6R0enVuZyBkZXMgTW9kZWxscy4g8J+YhAoKYGBge3J9CmZpdCA8LSBsbShlaW5rb21tZW4gfiBhbHRlciwgZGF0YSA9IGR0X3JhdykKdGlkeV9sbShmaXQpCmdsYW5jZV9sbShmaXQpCnBsb3RfbG0oZml0KQpgYGAKCgojIyBBbm5haG1lbiAoMS03KQoKMS4ga29ycmVrdGUgU3BlemlmaWthdGlvbiAoTGluZWFyaXTDpHQpCjIuIEVyd2FydHVuZ3N3ZXJ0IGRlciBSZXNpZHVlbiAobG9rYWwpID0gMAozLiBOb3JtYWx2ZXJ0ZWlsdGUgUmVzaWR1ZW4KNC4gVmFyaWFuemhvbW9nZW5pdMOkdCAoZ2xlaWNobcOkw59pZ2UgU3RyZXV1bmcgUmVzaWR1ZW4pCjUuIEJlcsO8Y2tzaWNodGlndW5nIGVpbmZsdXNzcmVpY2hlciBGw6RsbGUgKGtlaW5lIEF1c3JlaXNzZXIpCjYuIF9rZWluZV8gTXVsdGlrb2xsaW5lYXJpdMOkdCAoZ2VyaW5nZSBLb3JyZWxhdGlvbiB1bmFiaMOkbmdpZ2VyIFZhcmlhYmxlbikKNy4gX2tlaW5lXyBBdXRva29ycmVsYXRpb24gZGVyIFJlc2lkdWVuIChpbnNiZXMuIFplaXRyZWloZW4pCgoKIyMgS29ycmVrdGUgU3BlemlmaWthdGlvbgoKIyMjIERyaXR0LVZhcmlhYmxlbgoKIVtdKGJpbGRlci9icmllZnRyYWVnZXItZnJhdWVuLmpwZykKCkJyaWVmdHLDpGdlciDinInvuI8gYXJiZWl0ZW4gb2Z0IG51ciBoYWxidGFncyDwn5WQIHVuZCB2ZXJkaWVuZW4g8J+SsCBkYW5uIGF1Y2ggbnVyIGRpZSBIw6RsZnRlLgoKYGBge3J9CmR0IDwtCiAgZHRfcmF3ICU+JSAKICBtdXRhdGUoZWlua29tbWVuID0gaWZfZWxzZShpZCAlJSAzID09IDAsIGVpbmtvbW1lbi8yICsgcm5vcm0obiwgc2QgPSA1MCksIGVpbmtvbW1lbiksCiAgICAgICAgIGhhbGJ0YWdzID0gaWZfZWxzZShpZCAlJSAzID09IDAsICJqYSIsICJuZWluIikgJT4lIGZjdF9yZWxldmVsKCJqYSIsIGFmdGVyID0gSW5mKSkKCmdncGxvdChkdCwgYWVzKHggPSBhbHRlciwgeSA9IGVpbmtvbW1lbiwgY29sb3VyID0gaGFsYnRhZ3MpKSArCiAgZ2VvbV9wb2ludChhbHBoYSA9IDAuMikgKwogIHRoZW1lX3N0YXRhKCkKYGBgCgpPaG5lIERyaXR0LVZhcmlhYmxlIPCfmJ8KCmBgYHtyfQpmaXQgPC0gbG0oZWlua29tbWVuIH4gYWx0ZXIsIGRhdGEgPSBkdCkKdGlkeV9sbShmaXQpCmdsYW5jZV9sbShmaXQpCnBsb3RfbG0oZml0LCBidz01MCkKYGBgCgpNaXQgRHJpdHQtVmFyaWFibGUg8J+ZggoKYGBge3J9CmZpdCA8LSBsbShlaW5rb21tZW4gfiBhbHRlciArIGhhbGJ0YWdzLCBkYXRhID0gZHQpCnRpZHlfbG0oZml0KQpnbGFuY2VfbG0oZml0KQpwbG90X2xtKGZpdCkKYGBgCgpOb2NoIGJlc3NlciBtaXQgRHJpdHQtVmFyaWFibGUgdW5kIEludGVyYWt0aW9ucy1FZmZla3Qg8J+YgwoKYGBge3J9CmZpdCA8LSBsbShlaW5rb21tZW4gfiBhbHRlcipoYWxidGFncywgZGF0YSA9IGR0KQp0aWR5X2xtKGZpdCkKZ2xhbmNlX2xtKGZpdCkKcGxvdF9sbShmaXQpCmBgYAoKIyMjIEZ1bmt0aW9uYWxlIEZvcm0KCiFbXShiaWxkZXIvYnVuZGVzcG9zdC5qcGcpCgpCZWkgZGVyIEJ1bmRlc3Bvc3Qg8J+TryB3dXJkZSBlaW5lIEJlemFobHVuZyDwn5KwIG5hY2ggSHVtYW5rYXBpdGFsIGVpbmdlZsO8aHJ0LiAKCkVzIHd1cmRlIGVya2FubnQsIGRhc3MgYmVpIGp1bmdlbiBCZWFtdGVuIPCfp5EgZGllIFByb2R1a3Rpdml0w6R0IHN0w6Rya2VyIHfDpGNoc3QgYWxzIGJlaSDDhGx0ZXJlbiDwn6eTLiBEYWhlciBlcmhhbHRlbiBzaWUgaMO2aGVyZSBMb2huenV3w6RjaHNlIPCfkrAuCgpgYGB7cn0KZHQgPC0KICBkdF9yYXcgJT4lIAogIG11dGF0ZShlaW5rb21tZW5fcG9seW5vbSA9IC0zMDAgKyA4MCphbHRlciAtIDAuNzI1KmFsdGVyXjIsCiAgICAgICAgIGVpbmtvbW1lbiA9IGVpbmtvbW1lbl9wb2x5bm9tICsgcm5vcm0obiwgc2QgPSBzZF9laW5rb21tZW4pKQoKZ2dwbG90KGR0LCBhZXMoeCA9IGFsdGVyLCB5ID0gZWlua29tbWVuKSkgKwogIGdlb21fcG9pbnQoYWxwaGEgPSAwLjIpICsKICAgIGdlb21faml0dGVyKGFlcyh5ID0gZWlua29tbWVuX3BvbHlub20pLAogICAgICAgICAgICAgICAgd2lkdGggPSAwLjEsIGFscGhhID0gMC41LAogICAgICAgICAgICAgICAgY29sb3VyID0gImRhcmtibHVlIiwgc2hhcGUgPSAiLiIpICsKICB0aGVtZV9zdGF0YSgpCmBgYApFaW4gbGluZWFyZXMgTW9kZWxsIGlzdCBmYWxzY2ggc3BlemlmaXppZXJ0LiDwn5ifIERpZSBSZWdyZXNzaW9ucy1EaWFnbm9zZSB6ZWlndCBkaWVzLiDimJ3vuI8KCmBgYHtyfQpmaXQgPC0gbG0oZWlua29tbWVuIH4gYWx0ZXIsIGRhdGEgPSBkdCkKdGlkeV9sbShmaXQpCmdsYW5jZV9sbShmaXQpCnBsb3RfbG0oZml0KQpgYGAKCktvcnJla3QgaXN0IGVpbiBNb2RlbGwgbWl0IGVpbmVtIFBvbHlub20gKGhpZXIgZWluZW0gcXVhZHJhdGlzY2hlbSBUZXJtKS4g8J+YgwoKYGBge3J9CmZpdCA8LSBsbShlaW5rb21tZW4gfiBhbHRlciArIEkoYWx0ZXJeMiksIGRhdGEgPSBkdCkKdGlkeV9sbShmaXQpCmdsYW5jZV9sbShmaXQpCnBsb3RfbG0oZml0KQpgYGAKCiggRGVuIFIgSGlud2VpcyBhdWYgTXVsdGktS29sbGluZWFyaXTDpHQgaWdub3JpZXJlbiB3aXIgZWluZmFjaC4gTmF0w7xybGljaCBnaWJ0IGVzIGVpbmVuIFp1c2FtbWVuaGFuZyB6d2lzY2hlbiBgYWx0ZXJgIHVuZCBkZW0gcXVhZHJpZXJ0ZW4gYGFsdGVyXjJgLiApCgoKIyMgVmFyaWFuemhvbW9nZW5pdMOkdAoKIVtdKGJpbGRlci9taW5pc3Rlci5qcGcpCgpCZWkgZGVyIFBvc3Qg8J+TryBoYWJlbiBob2hlIEJlYW10ZSB1bmQgTWluc3RlciBlaW4gaMO2aGVyZXMgRWlua29tbWVuIPCfkrAgYWxzIGVpbmZhY2hlIEJlYW10ZS4gQWJlciBuaWNodCBhbGxlIHdlcmRlbiBob2hlIEJlYW10ZSDwn6eR4oCN8J+SvCBvZGVyIFtQb3N0LU1pbmlzdGVyXShodHRwczovL2RlLndpa2lwZWRpYS5vcmcvd2lraS9DaHJpc3RpYW5fU2Nod2Fyei1TY2hpbGxpbmcpIPCfj6QuCgpgYGB7cn0KZHQgPC0gZHRfcmF3ICU+JSBtdXRhdGUoZWlua29tbWVuID0gZWlua29tbWVuX2xpbmVhciArIGlkLzUwMCAqIHJub3JtKG4sIHNkID0gc2RfZWlua29tbWVuKSkKCmdncGxvdChkdCwgYWVzKHggPSBhbHRlciwgeSA9IGVpbmtvbW1lbikpICsKICBnZW9tX3BvaW50KGFscGhhID0gMC4yKSArCiAgICBnZW9tX2ppdHRlcihhZXMoeSA9IGVpbmtvbW1lbl9saW5lYXIpLAogICAgICAgICAgICAgICAgd2lkdGggPSAwLjEsIGFscGhhID0gMC41LAogICAgICAgICAgICAgICAgY29sb3VyID0gImRhcmtncmVlbiIsIHNoYXBlID0gIi4iKSArCiAgdGhlbWVfc3RhdGEoKQpgYGAKCkltIFN0YW5kYXJkLU1vZGVsbCBlcmhhbHRlbiB3aXIgdmVyemVycnRlIFN0YW5kYXJkLUZlaGxlci4g8J+YnwoKYGBge3J9CmZpdCA8LSBsbShlaW5rb21tZW4gfiBhbHRlciwgZGF0YSA9IGR0KQp0aWR5X2xtKGZpdCkKZ2xhbmNlX2xtKGZpdCkKcGxvdF9sbShmaXQpCmBgYAoKQmVzc2VyIGlzdCBlaW4gTW9kZWxsIG1pdCBfcm9idXN0ZW5fIPCfkqogU3RhbmRhcmQtRmVobGVybi4g8J+YhAoKYGBge3J9CmZpdCA8LSBsbV9yb2J1c3QoZWlua29tbWVuIH4gYWx0ZXIsIGRhdGEgPSBkdCkKdGlkeV9sbShmaXQpCmdsYW5jZShmaXQpICU+JSBtdXRhdGVfaWYoaXMubnVtZXJpYywgcm91bmQsIDIpIApgYGAKCiMjIE11bHRpa29sbGluZWFyaXTDpHQKCiFbXShiaWxkZXIvYnJpZWZ0cmFlZ2VyLWdydXBwZS5qcGcpCgrDhGx0ZXJlIPCfkbTwn5G1IFBvc3RiZWFtdGUg8J+TryBzaW5kIGtsZWluZXIgYWxzIErDvG5nZXJlIPCfkafwn5GmLgoKYGBge3J9CmR0IDwtIGR0X3JhdyAlPiUgbXV0YXRlKGdyb2Vzc2VfY20gPSByb3VuZChzZXEoMTgwLCAxNjAsIGxlbmd0aC5vdXQgPSBuKSArIHJub3JtKG4sIHNkID0gMi41KSwgMykpCgpnZ3Bsb3QoZHQsIGFlcyh4ID0gYWx0ZXIsIHkgPSBncm9lc3NlX2NtKSkgKwogIGdlb21fcG9pbnQoYWxwaGEgPSAwLjIpICsKICBnZW9tX3Ntb290aChtZXRob2QgPSBsbSkgKwogIHRoZW1lX3N0YXRhKCkKYGBgCkVzIGdpYnQgZGFubiBhdWNoIGVpbmVuIFp1c2FtbWVuaGFuZyDwn5OJIHp3aXNjaGVuIGBncm9lc3NlX2NtYCDwn5G1IHVuZCBgZWlua29tbWVuYCDwn5KwLCBuYXTDvHJsaWNoIGtlaW5lbiBrYXVzYWxlbiDimpnvuI8uCgpgYGB7cn0KZ2dwbG90KGR0LCBhZXMoeCA9IGdyb2Vzc2VfY20sIHkgPSBlaW5rb21tZW4pKSArCiAgZ2VvbV9wb2ludChhbHBoYSA9IDAuMikgKwogIGdlb21fc21vb3RoKG1ldGhvZCA9IGxtKSArCiAgdGhlbWVfc3RhdGEoKQpgYGAKCmBgYHtyfQpmaXQgPC0gbG0oZWlua29tbWVuIH4gZ3JvZXNzZV9jbSwgZGF0YSA9IGR0KQp0aWR5X2xtKGZpdCkKZ2xhbmNlX2xtKGZpdCkKYGBgCgpOZWhtZW4gd2lyIGBhbHRlcmAg8J+RtSB1bmQgYGdyb2Vzc2VfY21gIPCfk5AgaW4gZWluZW0gTW9kZWxsIPCfk4kgYXVmIGhhYmVuIHdpciBNdWx0aWtvbGxpbmVhcml0w6R0IPCflJcg8J+YnwoKYGBge3J9CmZpdCA8LSBsbShlaW5rb21tZW4gfiBhbHRlciArIGdyb2Vzc2VfY20sIGRhdGEgPSBkdCkKdGlkeV9sbShmaXQpCmdsYW5jZV9sbShmaXQpCnBsb3RfbG0oZml0KQpgYGAKCk11bHRpa29sbGluZWFyaXTDpHQgaW0gTW9kZWxsIPCfk4kgIGVya2VubmVuIPCflbXvuI8gd2lyIGF1Y2ggIGFuIGRlciBzZWhyIHN0YXJrZW4gS29ycmVsbGF0aW9uIGRlciB1bmFiaMOkbmdpZ2VuIFZhcmlhYmxlbi4KCmBgYHtyfQpkdCAlPiUgCiAgc2VsZWN0KGFsdGVyLCBlaW5rb21tZW4pICU+JSAKICBjb3IoKSAlPiUgCiAgYXNfdGliYmxlKCkgJT4lIAogIG11dGF0ZV9hbGwocm91bmQsIDIpCmBgYAoKCiMjIEF1c3JlaXNzZXIKCiFbXShiaWxkZXIvZmVybm1lbGRlYW10LmpwZykKCkJlaSBkZXIgRGF0ZW7DvGJlcnRyYWd1bmcg8J+ToSBnYWIgZXMgZWluZW4gRmVobGVyIOKaoO+4jyB1bmQgZGllIGxldHp0ZSBaaWZmZXIg8J+UoiB3dXJkZSBiZWkgZWluaWdlbiBHZWhhbHRzYW5nYWJlbiDwn5KwIGVudGZlcm50LiAKCmBgYHtyfQppZF9yYW5kb20gPC0gcnVuaWYoODAsIG1pbiA9IDEsIG1heCA9IG4pICU+JSByb3VuZCgpCgpkdCA8LQogIGR0X3JhdyAlPiUgCiAgbXV0YXRlKGVpbmtvbW1lbiA9IGlmX2Vsc2UoaWQgJWluJSBpZF9yYW5kb20sIHJvdW5kKGVpbmtvbW1lbi8xMCksIGVpbmtvbW1lbikpCgpnZ3Bsb3QoZHQsIGFlcyh4ID0gYWx0ZXIsIHkgPSBlaW5rb21tZW4pKSArCiAgZ2VvbV9wb2ludChhbHBoYSA9IDAuMikgKwogIHRoZW1lX3N0YXRhKCkKYGBgCgpXaXIgZXJoYWx0ZW4gZWluZSB2ZXJ6ZXJydGUgU2Now6R0enVuZy4g8J+YnwoKYGBge3J9CmZpdCA8LSBsbShlaW5rb21tZW4gfiBhbHRlciwgZGF0YSA9IGR0KQp0aWR5X2xtKGZpdCkKZ2xhbmNlX2xtKGZpdCkKcGxvdF9sbShmaXQsIGJ3PTUwKQpgYGAKCkRpZSBmYWxzY2hlbiBCZW9iYWNodHVuZ2VuIG3DvHNzZW4gd2lyIGVudGZlcm5lbiAob2RlciBrb3JyaWdpZXJlbikuIPCfmIMKCmBgYHtyfQpkdF90bXAgPC0gZHQgJT4lIGZpbHRlcihlaW5rb21tZW4gPiA1MDApCgpmaXQgPC0gbG0oZWlua29tbWVuIH4gYWx0ZXIsIGRhdGEgPSBkdF90bXApCnRpZHlfbG0oZml0KQpnbGFuY2VfbG0oZml0KQpwbG90X2xtKGZpdCkKYGBgCgoKCiMjIEF1dG9rb3JyZWxhdGlvbgoKIVtdKGJpbGRlci90ZWxlZm9uemVudHJhbGUuanBnKQoKRGllIFBvc3RiZWFtdGVuIPCfk68gd3VyZGVuIDIwMTAgdW5kIDIwMjAg8J+ThSBuYWNoIGlocmVtIEVpbmtvbW1lbiDwn5KwIGJlZnJhZ3QuIOKYju+4jwoKYGBge3J9CmR0X3RtcCA8LQogIGR0X3JhdyAlPiUgCiAgbXV0YXRlKGVpbmtvbW1lbiA9IHJvdW5kKGVpbmtvbW1lbiAtIHBhcmFtX2FsdGVyKjEwICsgcm5vcm0obiwgc2QgPSAxMCkpLAogICAgICAgICBhbHRlciA9IGFsdGVyIC0gMTAsCiAgICAgICAgIGJlZnJhZ3VuZyA9IDIwMTApICU+JSAKICBmaWx0ZXIoYWx0ZXIgPj0gMjApCgpkdCA8LQogIGR0X3JhdyAlPiUgCiAgbXV0YXRlKGJlZnJhZ3VuZyA9IDIwMjApICU+JSAKICByYmluZChkdF90bXApICU+JSAKICBhcnJhbmdlKGlkLCBiZWZyYWd1bmcpICU+JSAKICByZWxvY2F0ZShiZWZyYWd1bmcsIC5hZnRlciA9IGlkKQoKRFQ6OmRhdGF0YWJsZShkdCAlPiUgZmlsdGVyKGlkID4gNTAwKSkgICMgRGF0ZW4gYWIgSUQgNTAwIHVtIFBhbmVsLVN0cnVrdHVyIHp1IHZlcmRldXRsaWNoZW4gCmBgYAoKUG9zdGJlYW10ZSDwn5OvIGRpZSAyMDEwICDwn5OFIMOkbHRlciBhbHMgNTAgd2FyZW4gc2luZCBpbnp3aXNjaGVuIGluIFBlbnNpb24g8J+RtSB1bmQgaGFiZW4gYW4gZGVyIDIwMjAgQmVmcmFndW5nIOKYju+4jyBuaWNodCBtZWhyIHRlaWxnZW5vbW1lbi4KCmBgYHtyfQpkdF90bXAgPC0gZHQgJT4lIG11dGF0ZShiZWZyYWd1bmcgPSBmYWN0b3IoYmVmcmFndW5nKSkKCmdncGxvdChkdF90bXAsIGFlcyh4ID0gYWx0ZXIsIHkgPSBlaW5rb21tZW4sIGNvbG91ciA9IGJlZnJhZ3VuZykpICsKICBnZW9tX3BvaW50KGFscGhhID0gMC4yKSArCiAgdGhlbWVfc3RhdGEoKQpgYGAKCkltIE1vZGVsbCAobiA9IDE3NTApIHNpbmQgQmVvYmFjaHR1bmdlbiAodW5kIFJlc2lkdWVuKSBuaWNodCBtZWhyIHVuYWJow6RuZ2lnLiDwn5ifIFdpciBoYWJlbiBqYSB2aWVsZSB6d2VpbWFsIGJlZnJhZ3QuIOKYju+4jwoKYGBge3J9CmZpdCA8LSBsbShlaW5rb21tZW4gfiBhbHRlciwgZGF0YSA9IGR0X3RtcCkKdGlkeV9sbShmaXQpCmdsYW5jZV9sbShmaXQpCnBsb3RfbG0oZml0KQpgYGAKCiMgQXVmZ2FiZQoKIVtdKGJpbGRlci9yYWRpb3F1aXouanBnKQoKRGlza3V0aWVyZW4gU2llIPCfmYfwn5mH8J+ZhyBkaWUgUmVncmVzc2lvbnMtQW5uYWhtZW4gIPCfk4ggIHVuZCB6ZWlnZW4gU2llIG9iIHVuZCB3aWUgVmVybGV0enVuZ2VuIPCfj6UgZGVyIEFubmFobWVuIGVya2FubnQg8J+Vte+4jyAgd2VyZGVuIGvDtm5uZW4uCgotLS0KCiMgUXVlbGxlbgoKKyA8aHR0cHM6Ly9jb21tb25zLndpa2ltZWRpYS5vcmcvd2lraS9GaWxlOlNvcGhpZUFuZGVyc29uVGFrZXRoZWZhaXJmYWNlb2ZXb21hbi5qcGc+CisgPGh0dHBzOi8vY29tbW9ucy53aWtpbWVkaWEub3JnL3dpa2kvRmlsZTpCdW5kZXNhcmNoaXZfQmlsZF8xODMtNjMxNDItMDAwMSxfQ290dGJ1cyxfQnJpZWZ0ciVDMyVBNGdlcmlubmVuLmpwZz4KKyA8aHR0cHM6Ly9jb21tb25zLndpa2ltZWRpYS5vcmcvd2lraS9GaWxlOkJ1bmRlc2FyY2hpdl9CaWxkXzE4My0xOTkwLTA5MjQtMDIwLF9HZXJhLF9Qcm90ZXN0X3Zvbl9Qb3N0LUdld2Vya3NjaGFmdGVybi5qcGc+CisgPGh0dHBzOi8vY29tbW9ucy53aWtpbWVkaWEub3JnL3dpa2kvRmlsZTpCdW5kZXNhcmNoaXZfQmlsZF8xODMtVTAyMDUtMDAxMCxfQmVybGluLF9UYWdfZGVzX1Bvc3QtX3VuZF9GZXJubWVsZGV3ZXNlbnMsX1ZvcmJlcmVpdHVuZy5qcGc+CisgPGh0dHBzOi8vY29tbW9ucy53aWtpbWVkaWEub3JnL3dpa2kvRmlsZTpCdW5kZXNhcmNoaXZfQmlsZF8xODMtNTkxMTEtMDAwMixfQmVybGluLF9Qb3N0YW10X05fNTgsX0JyaWVmdHIlQzMlQTRnZXIuanBnPgorIDxodHRwczovL2NvbW1vbnMud2lraW1lZGlhLm9yZy93aWtpL0ZpbGU6VGVsZWZvbnplbnRyYWxlX2JlaV9kZXJfMi5JbmZhbnRlcmllZGl2aXNpb25fKEJpbGRJRF8xNTcwNTEyNSkuanBnPgorIDxodHRwczovL2NvbW1vbnMud2lraW1lZGlhLm9yZy93aWtpL0ZpbGU6UmFkaW9xdWl6XyUyMkFsbGVpbl9nZWdlbl9hbGxlJTIyX2F1c19kZW1fUmF0c3NhYWxfKEtpZWxfNDUuNTY3KS5qcGc+CgotLS0KCiFbXShiaWxkZXIvZmFpci5qcGcp